# This file esimates the effect conditional on cohabitation using only municipalities with recurring data

#  install.packages("xtable")
#  install.packages("ggplot2")
#  install.packages("rmeta")
#  install.packages("dplyr")
#  install.packages("rdrobust")

library(xtable)
library(ggplot2)
library(dplyr)
library(rmeta)
library(rdrobust)

## Define treatment dummy and assign weights 
load("het_data_recur.rdata")
data_both <- 
  data_coh_recur %>%
  mutate(treated = week > 0) 

## Loop over bandwidths

plot_data <- 
  data_frame(type  = rep(c("Conventional RD", "Robust RD"), each = 2),
             cohab = rep(c(0, 1), 2),
             eff   = rep(NA, 4),
             upp   = rep(NA, 4),
             low   = rep(NA, 4))

data_close <- 
  data_both %>%
  mutate(voted     = round(par_vote*n),
         abstained = round((1-par_vote)*n)) 

data_voted <- 
  as.data.frame(lapply(data_close, 
                       function(x,p) rep(x,p), 
                       data_close[["voted"]])) %>%
  mutate(voted = 1)

data_abstained <- 
  as.data.frame(lapply(data_close, 
                       function(x,p) rep(x,p), 
                       data_close[["abstained"]])) %>%
  mutate(voted = 0)

data_both <- 
  rbind(data_voted, data_abstained)

for(j in 1:2){
  data <- 
    data_both %>%
    filter(par_child_cohabit == (j-1)) 
  
  rd_2009 <- with(data %>% filter(year == 2009),
                  rdrobust(y = voted, x = week))
  sub_2009  <- c(rd_2009$coef[1], rd_2009$se[1])
  sub_2009_robust <- c(rd_2009$coef[3], rd_2009$se[3])
  
  rd_2013 <- with(data %>% filter(year == 2013),
                  rdrobust(y = voted, x = week))
  sub_2013 <- c(rd_2013$coef[1], rd_2013$se[1])
  sub_2013_robust <- c(rd_2013$coef[3], rd_2013$se[3])
  
  rd_2014 <- with(data %>% filter(year == 2014),
                  rdrobust(y = voted, x = week))
  sub_2014 <- c(rd_2014$coef[1], rd_2014$se[1])
  sub_2014_robust <- c(rd_2014$coef[3], rd_2014$se[3])
  
  rd_2015 <- with(data %>% filter(year == 2015),
                  rdrobust(y = voted, x = week))
  sub_2015 <- c(rd_2015$coef[1], rd_2015$se[1])
  sub_2015_robust <- c(rd_2015$coef[3], rd_2015$se[3])
  
  plot_data[j, 3] <- 
    meta.summaries(c(sub_2009[1],
                     sub_2013[1],
                     sub_2014[1],
                     sub_2015[1]),
                   c(sub_2009[2],
                     sub_2013[2],
                     sub_2014[2],
                     sub_2015[2]))$summary*100
  
  # store meta estimates. Point estimate plus upper and lower bound on 95% CI
  plot_data[j, 4] <- 
    plot_data[j, 3] + 
    meta.summaries(c(sub_2009[1],
                     sub_2013[1],
                     sub_2014[1],
                     sub_2015[1]),
                   c(sub_2009[2],
                     sub_2013[2],
                     sub_2014[2],
                     sub_2015[2]))$se.summary*100*qnorm(0.975)
  
  plot_data[j, 5] <- 
    plot_data[j, 3] + 
    meta.summaries(c(sub_2009[1],
                     sub_2013[1],
                     sub_2014[1],
                     sub_2015[1]),
                   c(sub_2009[2],
                     sub_2013[2],
                     sub_2014[2],
                     sub_2015[2]))$se.summary*100*qnorm(0.025)
  
  # Store robust estimates 
  plot_data[j + 2, 3] <- 
    meta.summaries(c(sub_2009_robust[1],
                     sub_2013_robust[1],
                     sub_2014_robust[1],
                     sub_2015_robust[1]),
                   c(sub_2009_robust[2],
                     sub_2013_robust[2],
                     sub_2014_robust[2],
                     sub_2015_robust[2]))$summary*100
  
  
  plot_data[j + 2, 4] <- 
    plot_data[j + 2, 3] + 
    meta.summaries(c(sub_2009_robust[1],
                     sub_2013_robust[1],
                     sub_2014_robust[1],
                     sub_2015_robust[1]),
                   c(sub_2009_robust[2],
                     sub_2013_robust[2],
                     sub_2014_robust[2],
                     sub_2015_robust[2]))$se.summary*100*qnorm(0.975)
  
  plot_data[j + 2, 5] <- 
    plot_data[j + 2, 3] + 
    meta.summaries(c(sub_2009_robust[1],
                     sub_2013_robust[1],
                     sub_2014_robust[1],
                     sub_2015_robust[1]),
                   c(sub_2009_robust[2],
                     sub_2013_robust[2],
                     sub_2014_robust[2],
                     sub_2015_robust[2]))$se.summary*100*qnorm(0.025)
}

plot_data <- 
  plot_data %>% 
  mutate(cohab = ifelse(cohab == 0, 0.5, 1),
         place = c(2.1, 1.1, 1.9, 0.9))

plot_cohab <- 
  ggplot(data = plot_data, 
         aes(x = eff, 
             y = place, 
             xmin = low,
             xmax = upp, 
             alpha = type)) + 
  geom_errorbarh(height = 0) + 
  geom_point(size = 2) + 
  scale_x_continuous("Effect on turnout in percentage points", 
                     breaks = seq(-6, 6, 2)) + 
  theme_classic() + 
  scale_alpha_discrete(name = "",
                       labels = c("Conventional RD", 
                                  "Robust RD")) + 
  scale_y_continuous("", breaks = c(1,2), 
                     labels = c("Parent does \nlive with child", 
                                "Parent does not \nlive with child")) +
  geom_vline(xintercept = 0, 
             alpha = 0.5, 
             linetype = "dashed") +
  theme(axis.text = element_text(size = 14),
        axis.title = element_text(size = 14),
        legend.text = element_text(size = 14),
        legend.title = element_text(size = 14))


